Chapter 5 Community composition

load("data/data.Rdata")

5.1 Taxonomy overview

5.1.1 Stacked barplot

genome_counts_filt_met<-genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  left_join(., sample_metadata, by = join_by(sample == EHI_number)) %>% #append sample metadata
  filter(count > 0) #filter 0 counts

genome_counts_filt_met$Elevation<-as.factor(genome_counts_filt_met$Elevation)
# Create an interaction variable for elevation and sample
genome_counts_filt_met$interaction_var <- interaction(genome_counts_filt_met$sample, genome_counts_filt_met$Elevation)

# Plot stacked barplot
ggplot(genome_counts_filt_met, aes(x=interaction_var,y=count,fill=phylum, group=phylum))+ #grouping enables keeping the same sorting of taxonomic units
    geom_bar(stat="identity", colour="white", linewidth=0.1)+ #plot stacked bars with white borders
    scale_fill_manual(values=phylum_colors) +
    labs(y = "Relative abundance", x="Elevation (m)") +
    guides(fill = guide_legend(ncol = 3)) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
          panel.background = element_blank(),
          panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black"),
          legend.position = "top",
    legend.title = element_blank(),
    legend.text = element_text(size=7))+
    scale_x_discrete(labels = function(x) gsub(".*\\.", "", x)) +
    facet_wrap(~Transect, scales = "free", labeller = as_labeller(function(label) gsub(".*\\.", "", label))) #only show elevation label

5.1.2 Phylum relative abundances

phylum_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>%
  left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
  left_join(genome_metadata, by = join_by(genome == genome)) %>%
  group_by(sample,phylum) %>%
  summarise(relabun=sum(count))

phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun),sd=sd(relabun)) %>%
    arrange(-mean) %>%
    tt()
tinytable_e6vorekz61jcz69pboiz
phylum mean sd
p__Bacillota_A 4.081110e-01 0.1600451966
p__Bacteroidota 3.813465e-01 0.1589155381
p__Pseudomonadota 5.849246e-02 0.0783694483
p__Bacillota 4.549049e-02 0.0635314708
p__Campylobacterota 3.481511e-02 0.0560037828
p__Desulfobacterota 2.666332e-02 0.0402945087
p__Verrucomicrobiota 1.518957e-02 0.0242135928
p__Bacillota_C 9.477405e-03 0.0104659743
p__Fusobacteriota 8.967841e-03 0.0368793328
p__Cyanobacteriota 3.965422e-03 0.0059815655
p__Spirochaetota 2.170591e-03 0.0096628708
p__Bacillota_B 2.069011e-03 0.0034571169
p__Actinomycetota 1.681425e-03 0.0051589130
p__Chlamydiota 1.027149e-03 0.0060969825
p__Deferribacterota 4.549927e-04 0.0026507460
p__Synergistota 7.766526e-05 0.0007996128 
phylum_arrange <- phylum_summary %>%
    group_by(phylum) %>%
    summarise(mean=mean(relabun)) %>%
    arrange(-mean) %>%
    select(phylum) %>%
    pull()

phylum_summary %>%
    filter(phylum %in% phylum_arrange) %>%
    mutate(phylum=factor(phylum,levels=rev(phylum_arrange))) %>%
    ggplot(aes(x=relabun, y=phylum, group=phylum, color=phylum)) +
        scale_color_manual(values=phylum_colors[rev(phylum_arrange)]) +
        geom_jitter(alpha=0.5) + 
        theme_minimal()

5.2 Taxonomy boxplot

5.2.1 Family

family_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == EHI_number)) %>% #append sample metadata
  left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,family) %>%
  summarise(relabun=sum(count))

family_summary %>%
    group_by(family) %>%
    summarise(mean=mean(relabun),sd=sd(relabun)) %>%
    arrange(-mean) %>%
    tt()
tinytable_6jgtb3jfy5kyj6p3lquy
family mean sd
f__Lachnospiraceae 2.825060e-01 0.1345643318
f__Bacteroidaceae 2.119504e-01 0.1049201320
f__Tannerellaceae 1.118797e-01 0.0753626317
f__ 6.502966e-02 0.0818702640
f__Helicobacteraceae 3.481511e-02 0.0560037828
f__Marinifilaceae 3.470570e-02 0.0255772494
f__UBA3700 2.788514e-02 0.0328903056
f__Desulfovibrionaceae 2.666332e-02 0.0402945087
f__Ruminococcaceae 2.425877e-02 0.0226024661
f__Rikenellaceae 1.885691e-02 0.0176408335
f__Erysipelotrichaceae 1.751136e-02 0.0268864962
f__Oscillospiraceae 1.683003e-02 0.0141001689
f__Coprobacillaceae 1.261022e-02 0.0335546952
f__Mycoplasmoidaceae 1.198842e-02 0.0269229159
f__Enterobacteriaceae 1.066950e-02 0.0648238061
f__Fusobacteriaceae 8.967841e-03 0.0368793328
f__Akkermansiaceae 8.595505e-03 0.0108643967
f__CAG-239 6.937215e-03 0.0097788795
f__LL51 6.594063e-03 0.0215416732
f__Anaerotignaceae 6.489923e-03 0.0073982451
f__UBA3830 5.877671e-03 0.0184934716
f__Gastranaerophilaceae 3.965422e-03 0.0059815655
f__Muribaculaceae 3.953836e-03 0.0407072381
f__Butyricicoccaceae 3.918292e-03 0.0050279924
f__CAG-274 3.022828e-03 0.0060610627
f__Acutalibacteraceae 2.697920e-03 0.0047318531
f__Anaerovoracaceae 2.658130e-03 0.0040447826
f__Pumilibacteraceae 2.520520e-03 0.0041575482
f__UBA1997 2.361296e-03 0.0079920448
f__CAG-508 2.302152e-03 0.0062849154
f__Brevinemataceae 2.170591e-03 0.0096628708
f__Peptococcaceae 2.069011e-03 0.0034571169
f__Rhodocyclaceae 1.918345e-03 0.0193815077
f__DTU072 1.775970e-03 0.0055679110
f__UBA660 1.726020e-03 0.0054429829
f__MGBC116941 1.698183e-03 0.0090561944
f__Massilibacillaceae 1.488714e-03 0.0024454858
f__Eggerthellaceae 1.364529e-03 0.0037021644
f__Enterococcaceae 8.375225e-04 0.0070514701
f__CALTSX01 5.149633e-04 0.0053018713
f__Chlamydiaceae 5.121857e-04 0.0030351423
f__CALVMC01 4.985785e-04 0.0019169750
f__Mucispirillaceae 4.549927e-04 0.0026507460
f__Clostridiaceae 4.173001e-04 0.0028996742
f__Acidaminococcaceae 3.966660e-04 0.0015706727
f__UBA1242 3.684261e-04 0.0014638710
f__RUG11792 3.682076e-04 0.0019335896
f__CAG-465 3.464214e-04 0.0015624302
f__Microbacteriaceae 3.168954e-04 0.0026872384
f__CAG-288 2.957044e-04 0.0017433720
f__Streptococcaceae 2.456082e-04 0.0018352415
f__Anaplasmataceae 2.412140e-04 0.0021407333
f__Hepatoplasmataceae 2.339936e-04 0.0024091117
f__Aeromonadaceae 2.305128e-04 0.0022380802
f__Peptostreptococcaceae 1.398982e-04 0.0014403401
f__Rhodobacteraceae 1.358917e-04 0.0013990908
f__Xanthomonadaceae 9.158653e-05 0.0009429410
f__Synergistaceae 7.766526e-05 0.0007996128
f__Lactobacillaceae 4.163513e-05 0.0004286599
f__Turicibacteraceae 0.000000e+00 0.0000000000

5.2.2 Genus

genus_summary <- genome_counts_filt %>%
  mutate_at(vars(-genome),~./sum(.)) %>% #apply TSS nornalisation
  pivot_longer(-genome, names_to = "sample", values_to = "count") %>% #reduce to minimum number of columns
  left_join(sample_metadata, by = join_by(sample == EHI_number)) %>% #append sample metadata
  left_join(genome_metadata, by = join_by(genome == genome)) %>% #append genome metadata
  group_by(sample,genus) %>%
  summarise(relabun=sum(count)) %>%
  filter(genus != "g__")

genus_summary %>%
    group_by(genus) %>%
    summarise(mean=mean(relabun),sd=sd(relabun)) %>%
    arrange(-mean) %>%
    tt()
tinytable_udj8m6zdn34yek68qhof
genus mean sd
g__Bacteroides 1.715631e-01 0.0995381723
g__Parabacteroides 1.099480e-01 0.0764703201
g__Roseburia 5.106281e-02 0.0740647328
g__Phocaeicola 3.848609e-02 0.0434625541
g__JAAYNV01 3.688241e-02 0.0750152647
g__Odoribacter 3.406510e-02 0.0255496567
g__Helicobacter_J 2.995383e-02 0.0398840003
g__CAG-95 1.693584e-02 0.0241924899
g__Alistipes 1.544141e-02 0.0156462125
g__Kineothrix 1.529781e-02 0.0582217172
g__MGBC136627 1.442882e-02 0.0229993736
g__Mycoplasmoides 1.128813e-02 0.0267998413
g__Hungatella_A 1.103483e-02 0.0667132423
g__Anaerotruncus 1.001128e-02 0.0112222704
g__Velocimicrobium 9.877812e-03 0.0159828933
g__Enterocloster 9.448113e-03 0.0098167225
g__Acetatifactor 9.119894e-03 0.0139402425
g__Fusobacterium_A 8.967841e-03 0.0368793328
g__Akkermansia 8.595505e-03 0.0108643967
g__Clostridium_Q 7.696833e-03 0.0137729617
g__Bilophila 7.309730e-03 0.0088486876
g__Lawsonia 7.063171e-03 0.0354503009
g__Intestinimonas 6.503115e-03 0.0061322665
g__Lacrimispora 6.449964e-03 0.0078252870
g__Lachnotalea 5.969286e-03 0.0091428683
g__Desulfovibrio 5.551727e-03 0.0084902145
g__MGBC140009 5.480884e-03 0.0133474239
g__Extibacter 5.251389e-03 0.0436025301
g__Coprobacillus 5.206125e-03 0.0159496122
g__Eisenbergiella 5.044262e-03 0.0077976628
g__Ventrimonas 5.035349e-03 0.0082923216
g__NHYM01 4.861275e-03 0.0425676661
g__Dielma 4.818065e-03 0.0065339333
g__CHH4-2 4.421332e-03 0.0044698946
g__RGIG4733 4.212224e-03 0.0102969273
g__Negativibacillus 4.095555e-03 0.0054792324
g__Thomasclavelia 3.821798e-03 0.0117039604
g__Hungatella 3.408253e-03 0.0046301753
g__C-19 3.303166e-03 0.0097049178
g__Citrobacter 3.169481e-03 0.0206845556
g__UMGS1251 2.856625e-03 0.0066227314
g__Oscillibacter 2.672493e-03 0.0038359772
g__CAZU01 2.662464e-03 0.0062677129
g__Copromonas 2.582309e-03 0.0037444602
g__Breznakia 2.545152e-03 0.0076365238
g__Mailhella 2.463379e-03 0.0034270857
g__Pseudoflavonifractor 2.327661e-03 0.0028208232
g__Intestinibacillus 2.304475e-03 0.0027604605
g__Escherichia 2.211570e-03 0.0159347159
g__MGBC165282 2.207931e-03 0.0052013921
g__Brevinema 2.170591e-03 0.0096628708
g__Rikenella 2.103581e-03 0.0033511480
g__Morganella 2.000896e-03 0.0206004803
g__Robinsoniella 1.982065e-03 0.0197231357
g__Parabacteroides_B 1.931687e-03 0.0064193758
g__Hafnia 1.921419e-03 0.0111820580
g__Fluviibacter 1.918345e-03 0.0193815077
g__JAIHAL01 1.891324e-03 0.0043767561
g__CAJLXD01 1.792508e-03 0.0041225034
g__Marseille-P3106 1.713488e-03 0.0025261999
g__UBA866 1.521215e-03 0.0026067620
g__MGBC116941 1.419756e-03 0.0090572569
g__Duncaniella 1.413153e-03 0.0145492957
g__RGIG6463 1.406224e-03 0.0032145982
g__Stoquefichus 1.405377e-03 0.0045901612
g__Limenecus 1.380444e-03 0.0029775228
g__JAAYQI01 1.259608e-03 0.0019959178
g__Lawsonibacter 1.166300e-03 0.0016508459
g__Scatousia 1.163255e-03 0.0033190366
g__MGBC101980 1.136457e-03 0.0043661490
g__Hespellia 1.092457e-03 0.0079658275
g__Clostridium_AQ 1.065062e-03 0.0036811111
g__Tidjanibacter 1.052202e-03 0.0029270036
g__Fournierella 1.012242e-03 0.0022501066
g__Eggerthella 9.852582e-04 0.0031664138
g__OM05-12 9.476666e-04 0.0022861996
g__CALXRO01 9.465351e-04 0.0060115481
g__CALURL01 9.003034e-04 0.0021568633
g__Harryflintia 8.821994e-04 0.0020357355
g__MGBC133411 8.788515e-04 0.0021894657
g__Scatacola_A 8.603252e-04 0.0028008223
g__Ventrisoma 8.513312e-04 0.0015698585
g__JALFVM01 8.313744e-04 0.0020032515
g__Bacteroides_G 8.279547e-04 0.0024893653
g__CAG-269 8.108577e-04 0.0029383065
g__IOR16 8.056128e-04 0.0024044902
g__CAG-873 7.720491e-04 0.0079487323
g__Buttiauxella 7.491009e-04 0.0061893836
g__14-2 7.174613e-04 0.0014179404
g__Ureaplasma 7.002851e-04 0.0022713191
g__Scatocola 6.837988e-04 0.0024038571
g__Dysosmobacter 6.815887e-04 0.0012774053
g__Muricomes 6.778722e-04 0.0023619925
g__Anaerovorax 6.716036e-04 0.0017450208
g__UBA7185 6.627870e-04 0.0018432854
g__Evtepia 6.474991e-04 0.0010541613
g__Butyricimonas 6.406058e-04 0.0016068308
g__MGBC131033 6.382237e-04 0.0015934066
g__CAJMNU01 6.222337e-04 0.0008887008
g__Beduini 5.858796e-04 0.0013098648
g__Muribaculum 5.497694e-04 0.0056602225
g__Scandinavium 5.341634e-04 0.0034686403
g__Lactonifactor 5.289757e-04 0.0013117146
g__CALTSX01 5.149633e-04 0.0053018713
g__CAG-485 5.079634e-04 0.0052298033
g__Merdicola 5.061702e-04 0.0018103952
g__Ventrenecus 4.907338e-04 0.0024331147
g__UMGS1202 4.839872e-04 0.0016896720
g__Copranaerobaculum 4.782162e-04 0.0029031775
g__NSJ-61 4.531589e-04 0.0013984436
g__Faecimonas 4.425060e-04 0.0017543086
g__RGIG8482 4.373523e-04 0.0020895489
g__Faecivivens 4.272845e-04 0.0008894505
g__RGIG9287 4.188304e-04 0.0020906022
g__Sarcina 4.173001e-04 0.0028996742
g__Blautia_A 4.105640e-04 0.0009994663
g__Scatenecus 4.064313e-04 0.0026411296
g__Phascolarctobacterium 3.966660e-04 0.0015706727
g__Raoultibacter 3.792710e-04 0.0011559737
g__Caccovivens 3.684261e-04 0.0014638710
g__CAJTFG01 3.655350e-04 0.0010180529
g__HGM11386 3.608509e-04 0.0015988134
g__CAG-465 3.464214e-04 0.0015624302
g__Amedibacillus 3.424477e-04 0.0023723526
g__Enterococcus_A 3.245757e-04 0.0023143225
g__UMGS2016 3.182000e-04 0.0012854921
g__Emergencia 3.175682e-04 0.0009661193
g__Holdemania 3.069052e-04 0.0010242364
g__Blautia 3.066198e-04 0.0011028540
g__Protoclostridium 3.038491e-04 0.0010789526
g__Fimivivens 3.014823e-04 0.0008070986
g__RGIG7389 2.957300e-04 0.0005828033
g__CAG-345 2.957044e-04 0.0017433720
g__UBA7173 2.935763e-04 0.0030225525
g__Bariatricus 2.900841e-04 0.0008344380
g__Agathobaculum 2.760992e-04 0.0016313981
g__CALXDZ01 2.635870e-04 0.0006108419
g__UBA940 2.597096e-04 0.0009035857
g__Microbacterium 2.538897e-04 0.0026139544
g__Aminipila 2.482408e-04 0.0007447481
g__Lactococcus 2.456082e-04 0.0018352415
g__Wolbachia 2.412140e-04 0.0021407333
g__Paramuribaculum 2.409489e-04 0.0024807206
g__Hepatoplasma 2.339936e-04 0.0024091117
g__Aeromonas 2.305128e-04 0.0022380802
g__WRHT01 2.165572e-04 0.0006925426
g__Zhenpiania 2.096592e-04 0.0012649339
g__UBA5026 2.092951e-04 0.0009734373
g__UMGS1663 1.980727e-04 0.0007253856
g__MGBC107952 1.736116e-04 0.0009841510
g__CALXEL01 1.709372e-04 0.0013664186
g__CAG-273 1.468577e-04 0.0007828475
g__Clostridioides 1.398982e-04 0.0014403401
g__Paracoccus 1.358917e-04 0.0013990908
g__JAAWBF01 1.274842e-04 0.0006116591
g__JAFLTL01 1.258715e-04 0.0012959262
g__Bacteroides_H 1.255258e-04 0.0012923672
g__RUG12867 9.167312e-05 0.0006100737
g__Stenotrophomonas 9.158653e-05 0.0009429410
g__Rahnella 8.286698e-05 0.0008021325
g__Lumbricidophila 6.300569e-05 0.0006486833
g__UBA3263 5.050541e-05 0.0005199850
g__Fructobacillus 4.163513e-05 0.0004286599
g__Clostridium 0.000000e+00 0.0000000000
g__Turicibacter 0.000000e+00 0.0000000000

5.3 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements (error in this section)
#Get list of present MAGs
present_MAGs <- genome_counts_filt %>%
    column_to_rownames(var = "genome") %>%
        filter(rowSums(.[, -1]) != 0) %>%
        rownames()

#Align KEGG annotations with present MAGs and remove all-zero and all-one traits
present_MAGs <- present_MAGs[present_MAGs %in% rownames(genome_gifts)]
genome_gifts_filt <- genome_gifts[present_MAGs,] %>%
            select_if(~!all(. == 0)) %>%  #remove all-zero modules
            select_if(~!all(. == 1)) #remove all-one modules

#Filter count table to only contain present MAGs after KEGG filtering
genome_counts_filt_filt <- genome_counts_filt %>%
  column_to_rownames(var = "genome")

genome_counts_filt_filt <- genome_counts_filt_filt[present_MAGs,]


dist <- genome_gifts_filt %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt_filt %>%
  #column_to_rownames(var = "genome") %>%
  #dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))
alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == EHI_number)) %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Transect, group=Transect, color=Transect, fill=Transect)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Transect",
          breaks=c("Aisa","Aran","Sentein","Tourmalet"),
          labels=c("Aisa","Aran","Sentein","Tourmalet"),
          values=c("#e5bd5b", "#6b7398","#e2815a", "#876b96")) +
      scale_fill_manual(name="Transect",
          breaks=c("Aisa","Aran","Sentein","Tourmalet"),
          labels=c("Aisa","Aran","Sentein","Tourmalet"),
          values=c("#e5bd5b50", "#6b739850","#e2815a50", "#876b9650")) +
      facet_wrap(. ~ metric, scales = "free", ncol=4) +
      coord_cartesian(xlim = c(1, NA)) +
      theme_classic() +
      theme(
        strip.background = element_blank(),
        panel.grid.minor.x = element_line(size = .1, color = "grey"),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)
      )

5.3.1 Regression plots

5.3.1.1 Richness diversity

columns <- c("richness","neutral","phylo","func","mapped","total")
alpha_div %>%
            select(sample,richness) %>%
            pivot_longer(-sample, names_to = "data", values_to = "value") %>%
            mutate(data = factor(data, levels = columns))   %>%
            left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
      ggplot(aes(x = Elevation, y = value)) +
        geom_point() +
        geom_smooth(method = lm) +
        facet_wrap(~ factor(Transect))+
        labs(x = "Elevation (m)")

5.3.1.2 Neutral diversity

alpha_div %>%
            select(sample,neutral) %>%
            pivot_longer(-sample, names_to = "data", values_to = "value") %>%
            mutate(data = factor(data, levels = columns))   %>%
            left_join(sample_metadata, by = join_by(sample == EHI_number))  %>%
      ggplot(aes(x = Elevation, y = value)) +
        geom_point() +
        geom_smooth(method = lm) +
        facet_wrap(~ factor(Transect))+
        labs(x = "Elevation (m)")

5.3.1.3 Phylogenetic diversity

alpha_div %>%
            select(sample,phylogenetic) %>%
            pivot_longer(-sample, names_to = "data", values_to = "value") %>%
            mutate(data = factor(data, levels = columns))   %>%
            left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
      ggplot(aes(x = Elevation, y = value)) +
        geom_point() +
        geom_smooth(method = lm) +
        facet_wrap(~ factor(Transect))+
        labs(x = "Elevation (m)")

5.3.1.4 Functional diversities

alpha_div %>%
            select(sample,functional) %>%
            pivot_longer(-sample, names_to = "data", values_to = "value") %>%
            mutate(data = factor(data, levels = columns))   %>%
            left_join(sample_metadata, by = join_by(sample == EHI_number))  %>%
      ggplot(aes(x = Elevation, y = value)) +
        geom_point() +
        geom_smooth(method = lm) +
        facet_wrap(~ factor(Transect))+
        labs(x = "Elevation (m)")

5.3.2 Mixed models

5.3.2.1 Richness diversity

alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
  lmerTest::lmer(richness ~ log(sequencing_depth) + Elevation*Transect + (1 | Sampling_point), data = ., REML = FALSE) %>%
  broom.mixed::tidy() %>%
  tt()

5.3.2.2 Neutral diversity

alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
  lmerTest::lmer(neutral ~ log(sequencing_depth) + Elevation*Transect + (1|Sampling_point), data = ., REML = FALSE) %>%
  broom.mixed::tidy() %>%
  tt()

5.3.2.3 Phylogenetic diversity

alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
  lmerTest::lmer(phylogenetic ~ log(sequencing_depth) + Elevation*Transect + (1 | Sampling_point), data = ., REML = FALSE) %>%
  broom.mixed::tidy() %>%
  tt()

5.3.2.4 Functional diversities

alpha_div %>%
  left_join(sample_metadata, by = join_by(sample == EHI_number)) %>%
  lmerTest::lmer(functional ~ log(sequencing_depth) + Elevation*Transect + (1 | Sampling_point), data = ., REML = FALSE) %>%
  broom.mixed::tidy() %>%
  tt()

5.4 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt_filt %>%
  #column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)